Q1

data <- readRDS("~/Downloads/week_11/data/historical_web_data_26112015.rds")
plot_ly(data, x = data$Latitude, y = data$Longitude, z = data$Depth, size = data$Magnitude) %>%  layout(scene = list(xaxis = list(title = 'Latitude'), yaxis = list(title = 'Longitude'), zaxis = list(title = 'Depth')))
## No trace type specified:
##   Based on info supplied, a 'scatter3d' trace seems appropriate.
##   Read more about this trace type -> https://plot.ly/r/reference/#scatter3d
## No scatter3d mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode

Q2

bigEQ = read_delim('~/Downloads/week_11/data/disaster.txt', '\t', escape_double = FALSE, trim_ws = TRUE)
sunami = bigEQ %>% filter(FLAG_TSUNAMI == 'Tsu')
sunami = sunami[!is.na(sunami$EQ_PRIMARY),]
sunami = sunami %>% filter(EQ_PRIMARY >= 8)
worldMap = map_data('world')
ggplot() + geom_polygon(data = worldMap, aes(x = long, y = lat, group = group),fill = 'blue') + geom_point(data = sunami, aes(x = LONGITUDE, y = LATITUDE), color = 'orange', size = 0.3) + transition_time(EQ_PRIMARY) + labs(title = 'mag {frame_time}') + ease_aes('linear')

Q3

iranEQ <- readRDS("~/Downloads/week_11/data/iran_earthquake.rds")
ggmap(read_rds("~/Downloads/week_11/data/Tehrn_map_6.rds")) + stat_density_2d(geom = 'polygon', data = iranEQ, aes(x = Long, y = Lat, fill = ..level.., alpha = ..level..))
## Warning: Removed 243 rows containing non-finite values (stat_density2d).

Q4

زلزله بزرگ را زلزله با قدرت بیشتر از ۶.۷ ریشتر در نظر گرفتیم. 

irBig = bigEQ %>% filter(EQ_PRIMARY > 6 & COUNTRY == 'IRAN')
distanceYear = data.frame(irBig[-1,]$YEAR - irBig[-nrow(irBig),]$YEAR)
colnames(distanceYear) = c('year')
notHappened = distanceYear %>% filter(year > 5)
cat(1 - nrow(notHappened)/nrow(distanceYear))
## 0.6764706

Q5

meannumeq <- bigEQ %>% group_by(COUNTRY) %>% summarise(numberOfEQ = n(), meanDeath = sum(DEATHS,na.rm = T)/numberOfEQ, alldeath = sum(DEATHS,na.rm = T))
mapDevice('x11')
heatMap = joinCountryData2Map(meannumeq, joinCode="NAME", nameJoinColumn="COUNTRY")
## 138 codes from your data successfully matched countries in the map
## 16 codes from your data failed to match with a country code in the map
## 105 codes from the map weren't represented in your data
mapCountryData(heatMap, nameColumnToPlot="alldeath")
## You asked for 7 quantiles, only 5 could be created in quantiles classification

Q6

fittedLine <- lm(bigEQ, formula = DEATHS ~  LONGITUDE + LATITUDE + FOCAL_DEPTH + EQ_PRIMARY )
summary(fittedLine)
## 
## Call:
## lm(formula = DEATHS ~ LONGITUDE + LATITUDE + FOCAL_DEPTH + EQ_PRIMARY, 
##     data = bigEQ)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -12556  -3519  -1872     35 311710 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -15451.303   3081.358  -5.014 6.13e-07 ***
## LONGITUDE       -1.232      5.813  -0.212  0.83226    
## LATITUDE        64.237     21.852   2.940  0.00335 ** 
## FOCAL_DEPTH    -21.929     11.354  -1.931  0.05367 .  
## EQ_PRIMARY    2678.740    464.175   5.771 1.01e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15440 on 1178 degrees of freedom
##   (4843 observations deleted due to missingness)
## Multiple R-squared:  0.031,  Adjusted R-squared:  0.0277 
## F-statistic:  9.42 on 4 and 1178 DF,  p-value: 1.7e-07

Q7

worldEQ = read.csv('~/Downloads/week_11/data/worldwide.csv')
worldEQ$date = as_date(worldEQ$time)
worldEQ$ID = interaction(day(worldEQ$date), year(worldEQ$date), month(worldEQ$date), sep='') 
worldEQ$country = sub(".*, ", "", worldEQ$place)
left = worldEQ %>% group_by(country, ID) %>% filter(mag == max(mag))
right = worldEQ %>% group_by(country, ID) %>% mutate(ind = rank(desc(mag))) %>% arrange(ind) %>% filter(ind == 2)
inter = inner_join(left, right, by = c('country', 'ID'))
train = inter[1:as.integer(0.9 * nrow(inter)),]
test = inter[-(1:as.integer(0.9 * nrow(inter))),]
model = lm(data = train, mag.x ~ mag.y)
summary(model)
## 
## Call:
## lm(formula = mag.x ~ mag.y, data = train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.39643 -0.28619 -0.09373  0.11327  2.61112 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.369494   0.022151   16.68   <2e-16 ***
## mag.y       1.005386   0.005515  182.31   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3763 on 7792 degrees of freedom
## Multiple R-squared:  0.8101, Adjusted R-squared:  0.8101 
## F-statistic: 3.324e+04 on 1 and 7792 DF,  p-value: < 2.2e-16
summary(stats::predict(model, test) - test$mag.x)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -2.2057 -0.1090  0.1338  0.0153  0.2921  0.3943

Q8

According to correlation coefficient and scatter plot, there is not any relation between depth and magnitude of earthquak

cor.test(worldEQ$depth,worldEQ$mag)
## 
##  Pearson's product-moment correlation
## 
## data:  worldEQ$depth and worldEQ$mag
## t = 50.03, df = 60629, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1914584 0.2067469
## sample estimates:
##       cor 
## 0.1991147

Q9

worldEQ = read.csv('~/Downloads/week_11/data/worldwide.csv')
worldEQ$country = sub(".*, ", "", worldEQ$place)
worldEQ$time = as.Date(worldEQ$time)
worldEQ$year = year(worldEQ$time)
worldEQ = worldEQ %>% group_by(place, year) %>% summarise(number = n()) %>% ungroup() %>% group_by(place) %>% summarise(n = mean(number))
worldEQ %>% arrange(desc(n)) %>% top_n(20) -> top20
## Selecting by n
ggplot(top20) + geom_bar(aes(x = place, y = n), stat = "identity")

Q10

bigEQ %>% arrange(desc(EQ_PRIMARY)) %>% top_n(10) ->f
## Selecting by TOTAL_HOUSES_DAMAGED_DESCRIPTION
ggplot(f) + geom_bar(aes(x = COUNTRY, y = EQ_PRIMARY),stat = "identity")
## Warning: Removed 3 rows containing missing values (position_stack).

bigEQ %>% filter(EQ_PRIMARY > 7) %>% group_by(COUNTRY) %>% summarise(n = n()) %>% top_n(10) -> d 
## Selecting by n
ggplot(d) + geom_bar(aes(x = COUNTRY, y = n),stat = "identity")

ts = sunami %>% filter(EQ_PRIMARY > 8)

ggplot() + geom_polygon(data = worldMap, aes(x = long, y = lat, group = group),fill = 'blue') + geom_point(data = ts, aes(x = LONGITUDE, y = LATITUDE), color = 'orange', size = 0.3)